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ABSTRACT 

We introduce a new approach to the modeUing of the hght distribution of galaxies, an orthonormal 
polar base formed by a combination of Chebyshev rational functions and Fourier polynomials that 
we call CHEF functions, or CHEFs. We have developed an orthonormalization process to apply this 
basis to pixelized images, and implemented the method as a Python pipeline. 

The new basis displays remarkable flexibility, being able to accurately fit all kinds of galaxy shapes, 
including irregulars, spirals, ellipticals, highly compact and highly elongated galaxies. It does this 
while using fewer components that similar methods, as shapelets, and without producing artifacts, 
due to the efficiency of the rational Chebyshev polynomials to fit quickly decaying functions like galaxy 
profiles. The method is lineal and very stable, and therefore capable of processing large numbers of 
galaxies in a fast and automated way. 

Due to the high quality of the fits in the central parts of the galaxies, and the efficiency of the 
CHEF basis modeling galaxy profiles up to very large distances, the method provides highly accurate 
estimates of total galaxy fiuxes and ellipticities. Future papers will explore in more detail the appli- 
cation of the method to perform multiband photometry, morphological classification and weak shear 
measurements. 

Subject headings: Methods: data analysis — Techniques: image processing — Galaxies: structure 
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1. INTRODUCTION 

Large surveys play an increasingly central role in Astronomy. The Sloan Digital Sky Survey (SDSS) (IMargonlll998D 
pioneered this trend, obtaining multi-color images covering more than a quarter of the sky, with a final dataset of 
230M celestial objects in an area of 8,400 square degrees, during its eight years of operations. The new generation 
of surveys is even more ambitious, and will generat e astonishing data volumes. The Panoramic Survey Telescope & 
Rapid Response System (Pan-STARRS) in Hawaii ([Kaiser et al.l 120021 ) . composed by four individual optical systems 
of 1.8 meter diameter each, will be able to survey 6,000 deg^ per night. It will catalog 99% of the stars in the 
northern hemisphere and make an ultradeep survey of 1200 square degrees to 17 = 27. The Dark Energy Survey (DES) 
(|Abbott et al.ll2005l ) will study the acceleration of the expanding universe by surveying 5,000 deg^ of the southern sky in 
five ban ds, and expects to det ect more than 300 million galaxies. The Javalambre-PAU Astrophysical Survey (J-PAS) 
project (jBemtez et al.|[2009aj) will survey 8,000 deg^ in 40 different bands, starting in 2013, measuring photometric 
redshifts with dz/{\ + z) ~ 0.003 and dz/{\-\ - z) ~ 0.01 for respectively, lOOM and 300M galaxies. The Large Synoptic 
burvey Telescope (LSST) (jlvezic et al.ll2008D wiU start to cover an area of 30,000 deg2 in 2015 with six optical fihers 
and expects to find a number of galaxies, 10 billion, dwarfing all other projects. To analyze and exploit such huge 
datasets, it is therefore necessary to develop techniques which are automated, objective and efficient. 

Modeling of the light distribution of galaxies is a crucial step for a wide array of problems in Galaxy Evolution 
and Observational Cosmology, which require accurate galaxy shape and flux measurements. The use of analytical 
proflles is the approach with longest tradition, employed e.g. by GIM2D (Simard 1998). This IRAF package analyzes 
galaxy images with a bulge/disk decomposition consisting of a Sersic profile (Sersic, .1968' ) plus an exponential, and 
optimization of the parameters is performed by the Me tropolis a l gorith m. Other works try to a djust galaxy profile s 
as a sum of elliptical Gaussians, as can be observed in IKuiikenI |l999l ) or in im2shape method ((iBridle et al.l l2003) . 
an d apply a maximum l ikelih ood method to fi t them . Lensing work have also contributed to develop new techniques, 
as llrwing fc Shmakoval (2005) or llrwing et al.l (|2007[ ). where a Gaussian times a polynomial profile is used as a model 
function. The idea of doing the model fitting in Fourier space has been also tried, for example in . Miller et al.l (|2007). 
being based on a de Vaucouleurs profile. 

The most widely used implementation of this approach is perhaps GALFIT (|Peng et al.ll2002l 120101 ) , which combines 
analytical profiles for the radial profile and trigonometric functions for the angular distribution. It works extremely 
well with objects with smooth features, as elliptical galaxies and thanks to its recently introduced angular components 
it displays considerable flexibility. However, as most analytical approaches, it requires considerable interactivity from 
the user in order to model complex galaxy proflles, especially at high S/N, since the number, type and initial values of 
the components have to be adjusted manually. Due to the large number of degrees of freedom, and the non-linearity 
of the fitting functions, which may require long convergence times or may not converge at all, it will be difficult to 
implement GALFIT, or any other analytical profile-based approaches, as a high-precision tool to handle the hundreds 
of millions of galaxies coming from the future large imaging surveys. 

A different approach, much better adapted to blind, automated fitting, is the use of orthonormal bases. The widely- 
known wavelets bases have been profusely applied to a great variety of science fields, including Astrophysics. They form 
a com plete and orthonormal system across scales, capable of a multiscale decompositi on of galaxies (Starck & Murtagfl 
|2002| ). They have inspired other bases, for example, beamlets (jDonoho fc Huoll2002| ) or curvelets fStarck et al..,2003( ). 
used to efficiently represent certain features present in astronomical images. 

The most succesful orthonormal basis methods in Astronomy are the so called shapelets (iMassey fc Refregieiil2005l : 
iRefregicr 2003). Buih using Hermite polynomials and Gaussians, these bases are extremely flexible, capable of decom- 
posing irregular objects and galaxies displaying small features. They are also fast and lineal. However, they feature an 
exponentially squared cut-off on the radial basis functions, which makes them fall significantly faster than most galaxy 
profiles. Because of this, they require ve ry large numbers of c oefficients to model objects with extended wings, what 
leads to oscillations in the radial shape (jMelchior et al.l[2010( ). This limitation seriously hampers their usefulness as 
automatic estimators of sh apes and fluxes for large galaxy surveys, a l though they have been used f or many astrophys- 
ical applications (jK uiikcn 1999; Rcfrcgier" 2003: Masscv et alj [200l iKellv fc McKavl 12004 120051: iGoldberg fc BaconI 
120051: iHevmans eTal. 2006: Massey et al. 2007; Kuiiken 2008(). An interesting attempt to tackle this issue are sersiclets 
( Ngan et al.l 120091) . a polar basis based on Sersic profiles. However, Sersic functions are poorly suited to form a 2D 
orthogonal basis, and this basis cannot fit irregular galaxy profiles. 

Given the obvious advantages of orthonormal bases, it seems highly desirable to develop one with the morphological 
flexibility and processing speed of the shapelets, but which is able to accurately and compactly represent realistic 
galaxy profiles up to very large radi i. Here we introduce such a basis, which we called CHEF functions or "CHEFs" , 
built using Chebyshev polynomials (|Bovdll 19*871) , which, as it is well known, disp lay remarkable propertie s of efficiency 
when fitting functions in the [-1,1] interval as their minimax property states (|Mason fc Handscombll2002[ ). The CHEF 
functions are separately built in polar coordinates, by expanding the radial component using Chebyshev rational 
functions (the composition of a Chebyshev polynomial and an algebraical mapping) and developing the azimuthal 
coordinate by means of sines and cosines. The result is a complete and orthonormal set of basis functions very well 
suited to fit astronomical objects able to fully model the whole flux of a galaxy to very large radii using a small number 
of coefficients. To explore the reliability and performance of this new technique, we have tested it on a wide sample 
of both simulated and real data, comparing the results with those obtained by the most popular publicly available 
techniques, GALFIT and IDL Shapelets software. 

The paper is organized as follows: in Sec. 2., we present the mathematical definition of the CHEF basis, demonstrate 
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the fast decay rate of the coefficients, what ensures the compactness of the basis, and derive the main galaxy parameters 
from the decomposition coefficients. Sec. 3. focuses on the practical implementation of the algorithm, describing how 
to discretize the basis and choose the free parameters of the fit, namely the scale size, the center and the decomposition 
order, and describes our functioning pipeline, which is being used to process data from different surveys. Sec. 4. shows 
the application of the method to different datasets, both real and simulated, and discusses its performance compared 
with that of shapelets. Sec. 5. states the main conclusions of the paper. 

2. MATHEMATICAL DEFINITIONS AND RESULTS 

We present an orthonormal basis set for image analysis based on polar basis functions that are constructed from 
rational Chebyshev functions, defined on the semi-infinite interval [0, +00), to expand the radial coordinate, and from 
sine and cosine waves to represent the azimuthal component. We call this basis "CHEF" functions (from H, the initial 
letter of Chebyshev, which is pronounced as "Che" in russian and F, the first letter of Fourier's family name) or, more 
compactly, "CHEFs". 

The Chebyshev polynomial of degree n is defined by the relation: 

Tn{z) — cos{n 9) , where 2; = cos 6*. (1) 

Chebyshev polynomials constitute a basis of the L^-space of squared-integrable functions defined on the finite 
interval [—1,1]. Using an algebraic coordinate transformation, this domain can be mapped onto the semi- infinite 
interval [0, -\rOo) and the rational Chebyshev functions can be obtained (|Bovdlll987f ): 



TL„(r) = T„ [y^j = cos \n ■ arccos [y^j j ■ (2) 

The constant L appearing in this expression will be called the scale parameter from now on and it is related to the 
"width" of the functions, and it does not affect their shape, as it can be seen in figure [TJ This scale is a free parameter 
of the decomposition and it has to be determined when applying the method to real data. Luckily, as we see below, 
the quality of the fit is not very sensitive to its exact value. 



Fig. 1. — Chebyshev rational function of order n = 5 with different scale sizes. 



These rational Chebyshev functions are orthogonal with respect to the weight function /(r) = \ — , that is: 

r + L y r 




TU{r-L)TLj{r-L)—-\-dr={Q, if i ^ j . (3) 



They tend monotously to 1 after their last minimum, more and more slowly as the order grows. Because of this 
property, and unlike shapelets, which fall off like a Gaussian, they can fit extended galaxy wings very compactly. At 
small radii, they quickly oscillate, and are thus capable of accurately describing the very fast change of fiux with 
radius displayed by most galaxies. The frontier between both regimes is roughly set by the scale parameter L. Their 
efficiency in modeling galaxy profiles is illustrated by figure [2] which shows the rational Chebyshev approximations to 
both a de Vaucouleur and a Sersic profile (with index = 0.75), using different Chebyshev rational functions orders. 
It can be observed how they accurately match both profiles over 3 orders of magnitude with less than 10 coefficients. 
In comparison, shapelets fail to model a de Vaucouleurs profile over two orders of magnitude even going to order 16. 
(see fig. 2 from (,Bosch.2010. )1. 
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Fig. 2. — Chebyshev approximations to both a de Vaucoulcurs (top panel) and a Scrsic profile (lower panel) (with index Us = 0.75), 
using Chebyshev rational functions of different order n. Order n < 10 is enough to provide high accuracy over the scales sampled by typical 
astronomical images. 



The polar CHEF functions are separable in r and 0, and built by expanding the radial component with normalized 
rational Chebyshev functions, whereas the angular coordinate is represented by a combination of sines and cosines: 

{^„,„(r,0;i)}„„ = TL^{r)Wmie)^ , (4) 

where C" = 2' if n > ' ^'^"^ T^m(^) is a general expression to represent both sin {mO) and cos {m9). 

CHEFs with scale size L = 1 and coefficient indices ranging from < n < 6 and —4 < to < 4 are displayed in figure 
[2 This set forms an orthonormal basis of the space of the finite-energy functions ([0, +00) x [—ir, tt] , (•, •)), which 
constitutes a Hilbert space with the inner product defined by 



+ 00 TT 

1 IL 



(Lg)^ J J fir,e) g{r,e)^^-dedr. (5) 
In this way, a smooth enough function / in polar coordinates can be decomposed into 

^ +OC +00 ^ +00 +00 

/ ^) = ^ E E ('^) ^™ (^) = ^ E E (™^) + fnn. TL^ (r) Sin (md)] (6) 

rn— n— m—0 n— 



and the coefficients fnm are calculated by 



TT 4-00 



fnm = I I fiz^4') TLn (z) ——j- \j - COS {m4>) dz d(j) 



27r2 7 J ' ' ' z + LS z 

TT +00 , (7) 

err 1 [l 



fnm ^ ^ J J fi^-'^) "^^'i (^) Y - sin {met)) dz d(l) 
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< n < 6 



Fig. 3. — Cosine and sine components of the first CHEF basis functions with scale size L = 1. 

These coefScients are finite, therefore, the series in 1^ is L^-convergent to the function / which is representing, 
what implies the expression in (|6]) is well defined. 

2.1. Coefficients decay 
CHEF coefficients are uniformly bounded by the norm of the function: 

\fnm\ = ^ |(/,0„™)| < ^ 11/11 • ||Ti„(r)W-„(0)|l = 11/11 (8) 
and their sum is also bounded 

m— n— 
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so the coefficients must tend to vanish as n and m increase. It can be proved (see appendix |Xl that the decay rate 
of these CHEF coefficients is 

|/«m| < -— — Vn,meN, (10) 
\n\ \m\ 2 

whenever the function / is smooth enough, and p is related to this smoothness; this ensures the quick convergence of 
the fits and therefore its compacteness. 

2.2. Shape parameters and fluxes 

The fact that CHEFs are radiafiy symmetric simpfifies the calculation of certain shape parameters as the total flux, 
the centroid, and the rms radius. As we will see, only the m G {—2, —1, 0, 1, 2} coefficients of this decomposition are 
included in the computation of these parameters. 

First of all, we define the integral = rL„(r) dr, with p>l and Re R. It can be proved that 
in the case n > 0, being 2F1 the hypergeometric function. If rt = 0, then simply 







P P+l 

Using this result, the fiux inside a circular aperture of radius R can be easily calculated from the CHEF cosine 
coefficients with to = 0: 

« - +00 



-ir 

In a similar way, the unweighted centroid {xc,yc) is 

R TT 

1 

Xc + tyc 



m) = / / fir, e)r dOdr ^2nJ2 fn,o ■ 



+00 

^ f [ f{r,e)r'{cose + ism9)dedr^^Y.[iJ^nS + fl-i)+^{fn,i-fn,-i)] I. 
»=o 



Finally, the expression for the quadrupole moments Fij = j j f{x,y)xiXjdxdy, yields this formula for the rms 
radius 



_ Fn + F22 _ 2^ 

rrms - p - p 

n=0 



— \ ' fc jn 

— p 2-^ Jn,Q ^3 ) 



and the ellipticity 



+00 

TT F _!_ 0,- F ^ (■^"^S + fn,2) ^3 

^11 — -^22 + ■i^i'l2 n=0 



Fii + F22 



Jn,0 ^3 

n=0 



Another definition of the ellipt icity, more useful as a shear estimator can be also obtained from the CHEF coefficients 
using the quadrupole moments (jMelchior et al.|[201Cl[ ): 

Fii ~ F22 + 2ii^i2 



Fii + F22 + 2(i^iii^22 ~ 


F^Y'^ 


+00 

J2 {fn.2 + fn.2) F^ 
ji=0 










+00 / +00 

2E/V3+1/ E 

Tl = V iiJ2=0 


{fl\,0 + -2ifL-2 + 




\(fh,2 f 


\-2)ifl2,2 f 


2,-2) 


^3 ^3 



(11) 
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Using a different approach, the behaviour of the CHEF coe fScients under linear transformations described in appendix 
[Cl the asymmetry parameter described by iConselicd (|2003D can be calculated using the CHEF coefficients; 



-oo +00 



E E + sin(m0)) TL„(r) r d^dr 

m— O,odd 71—0 



R2 

And, finally, the Sersic index can be also easily expressed by means of the CHEF coefficients: 

— fc In (r/rp) 

" = 7 Z — 1 -^-^ 

/ +00 -\-oo 

In E E [/^™cos(™0) + /^„sin(m0)] TL„(r) 

\ n— m— 

where is the effective radius of the galaxy, Se the surface brightness at r^, and fc is a constant coupled to n such 
that half of the total flux is within r^- 

3. PRACTICAL IMPLEMENTATION 

In this section we also show how to deal with the free parameters of the decomposition, namely the scale parameter 
L, the center (xc, j/c) and the total number of coefficients N and Af , something necessary for any practical application. 

3.1. Discrete basis functions 

Astronomical images are formed by pixels, and therefore, any analytical basis must be discretized to model real 
data, going from the space of continuous functions defined on to the space of matrices with size pi 'x P2, that is, 
the size of the image being fit. This is achieved by integrating the CHEF basis within each pixel: 

{0„™(r„0fc,L)}_- <[ JJ ^TL^{r,)Wrn(Ok)} , (12) 

(^(rj.Sfc) pixel 

where [vj, 0k) are the coordinates of the nodes of a polar grid of size pi x p2. 

It is obvious that this step will not produce a complete, let alone orthonormal, basis on the pixel space. In principle, 
one could sample the basis at the pixel centers and fit the galaxies using an optimization scheme. The basis is so 
well suited to fit galaxy profiles that it can be used in this way with reasonably good results. However, completeness 
and orthonormality are essential for a robust and automated approach. The first property ensures the flexibility of 
the basis, the capability of fitting all possible galaxy profiles, including irregulars. The second is crucial for a fast 
implementation, since the coefficients of the decomposition can be calculated in a straightforward way, using the inner 
product defined in the space of real vectors. 

Therefore, to orthonormalize our discrete basis, we chose the Modified Gram-Sch midt method (MGS). whic h consists 
in a smart rearrangement of the calculations of the Classical Gram-Schmidt process (IGolub fc Van Loanl[l996h . yielding 
the same results as CGS in exact arithmetic, but displaying a much more stable behavior in finite-precision arithmetic. 
In fact, due to rounding errors, the application of the classic method produces final vectors which are o ften not 
orthogonal. A description of the MGS computation algorithm can be looked up in (|Golub fc Van LoanlHoOfih . When 
applied to our pixelized bases, it produces nicely orthonormal vectors, up to the larger orders tested. Using the new 
basis (pnm we have that the image light distribution / can be decomposed as 

NM 

f{r],Sk) = ^ fnm(f>nrn{rj,9k). (13) 
nrn 

The values of the decomposition coefficients are fast and accurately calculated by means of the usual inner product 
in the matrix space, that is: 

fnm = '^f{rj,ek)(f>nrn{rj,Ok), (14) 

where the sum j, goes over all the pixels in the image. 

The orthonormalization method needs to be applied only once, for a given scale and image size. Given that the CHEF 
reconstruction is relatively robust to the scale L and image size, it is possible to produce a grid of basis functions, with 
different scales and sizes and precalculate the orthonormalization in advance. This significantly improves the efficiency 
of the method when applied to large data volumes. 

3.2. Number of coefficients, center and scale size 

As in the case of any similar method, the decomposition has several free parameters: the scale L, the total number 
of coefScients (A^ -I- 1) • (2M -|- 1) and the origin {xc, TJc) ■ It is obvious that if we grossly err in the choice of one of them. 
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the algorithm wih not work properly. Conveniently, quasi-optimal choices for all these parameters, except N and M, 
can be obtained directly from the output of SExtractor. 

The center of the reconstruction is a sensitive parameter, since the basis functions have an extreme at its position: 
it should match the maximum of the light distribution of the object, since a wrongly determined center would require 
an innecessarily high number of coefficients for an accurate fit. Luckly the centroid calculated by SExtractor provides 
a suitable center in practically all cases. 

The scale parameter L is not related to the shape of the CHEF basis, but to their width (the scale has an "accordion" 
effect over the basis functions). The larger the scale is, the slower the basis functions reach their last extreme, as can 
be seen in figure [TJ A scale size which is too small would overfit the high frequency noise in the central part of the 
image, whereas a scale size which is too large would smooth over the galaxies, incapable of fitting the central peak 
and the tiny features. In principle one could perform the fits iteratively, looking for the best scale L. However, we 
have verified that the quality of the fit is robust with respect to the value of L and using the SExtractor half-light 
radius estimate as a proxy for its value yields extremely accurate results. It is therefore not necessary to optimize this 
parameter. 

The remaining parameter, the orders of the decomposition N and M need to be carefully selected. We must choose 
a number of coefficients which is large enough to accurately model the image, but without overfitting the noise. The 
number of coefficients is thus optimized until the reconstruction residual approaches 1: 



where / is the observed image, fnm are its CHEF coefficients, is referred to the pixel noise map, Upi^ei is the 
number of pixels in / (that is, ripixei — Pi ■ P2), and ricoeff represents the number of CHEF coefficients used to 
decompose / (i.e., ricoeff — {N + 1){2M + 1)). At the point where ~ Ij the residual between the observed data 
and the model is consistent with the noise level. 

3.3. The CHEF pipeline 

The process to decompose real image data using the C HEF basis has been implemented in a P ython software package, 
which is being applied to the ALHAMBRA survey data (jMoles et al.ll2008l : lBenftez et al.ll2009b[ ). In the process outlined 
below we do not deal separately with the PSFs, fitting directly the objects without attempting any deconvolution. 
That will be the subject of a separate paper. 

1. Image analysis with SExtractor and source filtering. We run SExtractor on the image, to select galaxies and 
estimate the parameters we need for the fit, namely the centroid and the half-light radius, which will be used to 
determine the origin and the scale size of the CHEFs, respectively. 

2. Choice of the reconstruction size and the maximum number of coefficients. We have found that a good choice 
for the size of the reconstruction is to set Lmax = 4 x rhi, where r^i is the half-light radius of the galaxy, as 
measured by SExtractor. In practically all cases this fully includes the galaxy flux. 

As a general rule, the larger the number of pixels in the image to be fit, the larger the required number of 
coefficients. For instance, for two images of the same galaxy, the one with a smaller pixel scale will require 
more coefficients (assuming of course that the S/N is similar in both cases). We set 7V""'^,M""'^ and choose 
this value according to the image size. For instance, for a regular image of 60 x 60 pixels the usual choice is 
js^rnax ^ ^jniax ^ 7^ whcrcas for morc extended sources of about 120 x 120 pixels we set N"'"'' = M"""^ = 10. 

Once the image size is determined, we mask out all the pixels belonging to other objects present within the 



3. Evaluation of the basis. The CHEF basis corresponding to N"^"-^ and M™°^ is evaluated and orthonormalized 
using the modified Gram-Schmidt algorithm (see section [OJ • Thanks to the robustness of the CHEF decompo- 
sition, it is possible to precalculate this step, for a grid of different scale parameters L, so that the software just 
has to look for the closest value stored in the database, making it computationally much faster. For large data 
volumes this results in a speed up of more than a factor of 3. 

4. Determination of the optimal number of coefficients. As explained in Sec. 13.21 we choose the optimal number of 
coefficients N and M iteratively, by choosing the values which yield ~ 1- The corresponding basis functions, 
up to the required order, are extracted from the whole orthonormalized set previously computed in step 3. 




(15) 



frame. 



5. 



Evaluation of CHEF coefficients and calculation of the model. The Ucoeff — [N + 1) x (2M + 1) coefficients of 
the CHEF expansion and the resulting model are calculated by simple scalar products, as described in equations 
(HSl) and (O. 
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(c) 

Fig. 4. — CHEF decomposition of a spiral galaxy from the UDF: a) original image, b) CHEF reconstruction using N 
M = 14 angular coefficients, (c) residuals 



15 radial and 



(a) 



(b) 



(c) 



Fig. 5. — CHEF decomposition of an edge-on galaxy from the UDF: a) original image b) CHEF reconstruction with N 
and c) residuals (note that the small object below the main galaxy is also fit separately) 



: 5 and M = 7, 



6. Iteration of the algorithm for overlapping objects. For those objects whose flux may be affected by nearby objects 
we repeat steps 2-5, but instead of masking out the objects, we subtract their profiles already calculated in the 
previous steps. 

Figures mSl and |6] show some examples of Hubble Ultra Deep Field (jBeckwith et al.l[2004t iCoe et al.l 12006 1) galax- 
ies. The large amount of detail revealed in those images makes them a useful benchmark for galaxy decomposition 
methods. These objects, while typical, are not trivial to fit due to the complicated structure revealed by the high S/N 
observations. The first is a spiral galaxy, displayed in figure^ The image frame has 121 x 121 pixels, and the algorithm 
establishes that the best fit is obtained with = 15 and M — 14, i.e., using ricoe// — (15 + 1)(14 x 2 -|- 1) = 464 
coefficients. The very low level of the residuals shows that the CHEF basis can easily fit both the central bulge and 
large structures like spiral arms. We apply the same procedure to other two galaxies, an edge on galaxy with very 
high observed ellipticity (Fig. [5]) and a highly irregular galaxy (Fig. [6]). The fits are also excellent, what illustrates 
the flexibility of the CHEF basis, able to handle most galaxy morphologies. 

4. COMPARISON WITH OTHER METHODS 

We compare the C HEF decomposition with two of the most widely available techniq ues for galaxy modeling, GAL- 
FIT (jPeng et al.l [2002f ) and the IDL Shapelets software. GALFIT (jPeng et al.l 12003) uses several two-dimensional 
analytic models to fit the galaxy images, as the Sersic profile, the exponential disk, the Nuker law or the Gaussian 
and Moffat /L orentzian fun ctions. It is possible to use an arbitrary number of components. The shapelets software 
assey fc Refregiei] 120051 ) is based on an orthonormal basis built from Hermite polynomials and Gaussians. Figure 
[7] shows the results of fitting 12 randomly chosen galaxies from UDF with GALFIT, IDL Shapelets and our CHEF 
pipeline. 

The GALFIT results show that radially symmetrical multicomponent models fail to adequately represent the wide 
variety of shapes displayed by galaxies observed at high S/N, and that it is essential to include angular components, 
as indeed the new version of GALFIT does. However, an essential problem of GALFIT, and in general of any similar 
approach, is that by using a combination of several analytical profiles, GALFIT is mathematically equivalent to a 
decomposition into a set of functions which not only do not form an orthogonal basis, but are highly degenerate 
among themselves. Therefore as soon as the number of components increases, the fits must be performed interactively, 
in order to provide the right initial parameters, the number and type of the components, etc. This makes very difficult 
its use to obtain automatic, highly accurate fits of large numbers of galaxies. 
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(a) 



(b) 




Fig. 6. — CHEF decomposition of an irregular galaxy from the UDF: a) original image b) CHEF reconstruction with N = 5 and M - 
and c) residuals 



We have run the IDL shapelets software allowing Umax = 20 as the maximum order, equivalent to Ucoeff = 400 
coefficients. For comparison, we use TV™"^ = f^/jmax _ ^ equivalent to ricoeff = 231 coefficients, for the CHEF 
decomposition, which also runs faster, even taking into account the orthonormalization process. Despite of using more 
coefficients, we see the shapelet reconstruction often fails to adequately fit all galaxy shapes and t ends to produce 
ring- like artifacts, probably as a consequence of the inestability of the radial proffie fits ()Boschll2010l ). 

The CHEF decomposition manages to keep the residuals to a very low level, fitting very complicated galaxy shapes, 
producing very few artifacts. In addition, as we will see below, it obtains highly accurate estimates of galaxy total 
fluxes and shapes. 
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Fig. 7. — Reconstructions and obtained residuals obtained after applying the three metliods to an ensemble of galaxies randomly extracted 
from UDF: images on the top correspond to the original galaxies, first row below it represents the reconstructions obtained by shapelets, 
GALFIT and CHEF softwares (from left to right), and second row shows the resulting residuals, in linear scale. 
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Fig. 8. — Ellipticity absolute error obtained by applying shapelet and CHEF algorithms to a sample of simulated galaxies with different 
Sersic profiles and ellipticity parameters. 
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Fig. 9. — Magnitude relative error obtained by applying shapelet and CHEF algorithms to a sample of simulated galaxies with different 
Sersic profiles and ellipticity parameters. 



4.1. Magnitudes and ellipticities 

As stated above, the main motivation to devefop the CHEF basis is having a tool to automatically measure accurate 
fluxes and shapes in the context of large galaxy surveys. 

Here we test the performance of CHEFs using an ensemble of 350 mock galaxy profiles, simulated using Sersic profiles 
with indices ranging from 0.5 to 4, and sheared by different amounts, ranging from ellipticity to 0.5, by applying to 
the cartesian coordinates the following transformation: 



1 - 71 -72 
-72 1 + 71 



We then fit the galaxy shapes using the IDL shapelet software and our CHEF pipeline. We can see in Fig |5| that, 
as expected, shapelets fail to adequately objects with high ellipticities, specially for Sersic indexes n w 1.2, where 
the relative error in the measurement of the ellipticity reaches 10%. Although we are using circular shapelets here, it 
seems that elliptical shapelets are also affected by this problem (Bosch 2010). 

Using the formula (ITT]) , we estimate the ellipticities of our CHEF fits. Our error is < 0.5% for all the explored com- 
binations of profiles and shears. This shows that CHEFs are very promising as a tool for weak lensing measurements, 
and issue which will be explored elsewhere. 

To verify the performance of the CHEFs for photometric purposes, we added Gaussian noise to the galaxies above, 
and applied both the shapelets and CHEF algorithms to the data. The results are presented in Fig [HI Again we see 
that CHEFs provide highly accurate measurements of the total flux for almost all combinations of proflle steepness 
and ellipticity, with a typical error of only « 0.4%. Shapelets, on the other hand, introduce up to 4% errors in this 
very high S/N simulated images. 
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Motivated by the need to perform accurate galaxy photometry and shape measurements in large surveys like the 
upcoming JPAS (|Bemtez et al.|[2009af) , we have developed a new orthonormal basis formed by a combination of rational 
Chebyshev and trigonometric functions for, respectively, the radial and angular components, and that we called CHEF 
functions, or CHEFs. 

The new basis displays remarkable flexibility, being able to accurately fit all kinds of galaxy shapes, including 
irregulars, spirals, ellipticals, highly compact and highly sheared galaxies. It does this while using fewer components 
that similar methods, as shapelets, and without producing artifacts, due to the efficiency of the rational Chebyshev 
polynomials to fit quickly decaying functions like galaxy profiles. The method is lineal as well as very stable and 
robust, therefore, able to work very fast and to automatically process large numbers of galaxies. 

Due to the high quality of the fits in the central parts of the galaxies, and the ability of the CHEF basis to compactly 
model galaxy profiles up to very large distances, the method provides highly accurate estimates of total galaxy fluxes 
and ellipticities. 

Future work will explore in more detail the application of the method to problems like multiband photometry, weak 
shear measurements, deconvolution, morphological classification, etc. 

The authors thank Renato Alvarez for reading an early draft and making useful comments. We also had interesting 
and useful conversations with Sarah Bridle, Peter Melchior and Joel Berge. We acknowledge the support and help 
received during the visits to Jet Propulsion Laboratory, University College of London and Space Telescope Science 
Institute. This work was made possible by the financial support of the Consejo Superior de Investigaciones Cientificas 
and its ISP grants program, the ALHAMBRA grant AYA2006-14056 and the government project AYA2010-22111- 
C03-00. 

APPENDIX 
COEFFICIENTS DECAY 

Let us calculate the decay rate of the CHEF coefficients using the fact that a Chebyshev series is a Fourier series 
with a change of variable, in such a way that the Chebyshev coefficients of a univariate function g are proportional to 
the Fourier coefficients of the composition g o cos: 

gn = 2G[n], Vn G N, with G(i) = .g(cos(i)). (Al) 
Proposition A.l. Let f be a function of the space L^{[-~1, 1]). Let us define the function F in the following way: 

(A2) 

t i-> f[cost) ^ ' 

Then, the Fourier coefficients of F are proportional to the Chebyshev coefficients of f , that is, 

i?W-i/„, VneN (A3) 

where fn is the n-th coefficient in the Chebyshev series of f and ^^j^^ — 2' i/n>0' 
Proof : The function is a 27r-periodic function, so its Fourier series can be calculated: 

7r 7r 

F[n]^— I F{t)e-"'^ dt ^ — ( i^(i)e-"* dt + — [ i^(i)e-*"* dt = 
2vr J 2t: J 2Tr J 

— 7r — 7r 



= — / F(t) re™*+e-™*l dt=- I F(t]cos(nt)dt 

271 J TT J 



1 1 



= - / f{x)cos{n-arccoax){l~x^)-^/^dx^- [ f{x)T^{x){l - x^y^^^ dx = y fn- 

TT J TT J k 



This result is extremely useful since it allows us to transfer the Fourier series properties to the CHEF one. In addition 
to this, it is well-known the Fourier coefficients of a smooth enough function / are bounded, as the following theorem 
states. 

Theorem A. 2. Wizel^VlOSA) Let f be a periodic function with an interval of periodicity [0,P), continuous with its 
derivatives up to order p — 1, p > 1, that is, f € C^~^ . Let its derivative of order p be piecewise continuous in the 
interval. Then, there exists a constant A > such that, for Fourier coefficients of the function f , we have 

|/[/c]|<^, VfceN. (A4) 
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Therefore, using (jA3l) and (|A4I) it is easy to reach the decay rate of the CHEF coefficients, as it is shown in the 
following proposition. 

Proposition A. 3. Let f be a function of the space £^([0, +oo] x [— tt, tt]) defined in polar coordinates. Let us assume 
that f is continuous in the angular coordinate with its partial derivatives up to order p — 2, and that its partial 
derivative of order p — 1 is piecewise continuous in the interval [— 7r,7r]. Let f be also piecewise continuous in the 
azimuthal coordinate (the interval [0, +oo)). Then, there exists a constant C > such that for the CHEF coefficients 
of f it is had the following boundary: 

C 

l/«m| < --— TTT, Vn,meN. (A5) 
\n\ \m\ 2 

Proof : Let us define the function g as 

g : [-7r,7r] M 



^ / f{r,e)TLr,{r)-^J-dr 
J r+L\r 





This function can be extended in a periodic way it fulfils that it is continuous with its derivatives up to order p — 1, 
as well as the p-th derivative is piecewise continuous in the periodicity interval tt, vr]. Then, using the theorem I A. 21 
there exists a non-negative constant A such that 

Moreover, there exists a relationship between the CHEF coefficients of / and these Fourier coefficients of g: 

OO TT , TT 

fnm ^ J J fir, 0) TL4r) e-™« -^^^ dr ^ J g{0) e-""« de = 27rg[m]. 

— TT — TT 

Therefore, 

On another front, let us define a function h as 

h: [-7r,7r] R 

TT 

t ^ I f e-™« de ■ 

Let us consider the function h — ho cos. This function can be extended in a periodic way to M and it is continuous 
within its inter val o f periodicity. Its derivative is piecewise continuous in this interval, so we are under the assumption 
of the theorem IA.2I and it can be stated that there exists a constant B > such that 



Working in a similar way, we reach the following relation 

1 



fn 



f ( L, 9 ] c-""« d9 



-1 

1 



T„(z)(l-z2)-i/2dz = /i„. 



Using the proposition lA.il we have that h[n] = — hn, then 

k 



\fnm\ < 4- (A7) 



Finally, using (|A6|) and (|A7p . it is obtained 



If P< '^'^ 

\Jnm\ _^ 



p+l„2 • 
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STURM-LIOUVILLE EQUATIONS 
The so called Sturm- Liouville differential equations fit the following expression 

[r{x)y']' + [qix)+Xp{x)]y^O. (Bl) 

The solutions y{x) of this ordinary differential equation are called eigenf unctions while the valid values A associated to 
them receive the name of eigenvalues. It can be proved that if the functions p, q and r are real and continuous over an 
interval [a, b] being p{x) > (or p{x) < 0),Va; G [a, b], then all the eigenvalues lambda are real and the eigenfunctions 
y are orthogonal in [a, b] with the weight function p. 



It is easy to see that the rational Chebyshev functions constitutes a family of eigenfunctions of a Sturm-Liouville 
equations, specifically the {T Ln{x)} ^ functions form the solutions set of the following EDO 

displaying, therefore, all the properties of the Sturm-Liouville eigenfunctions. 



LINEAR TRANSFORMATIONS 

It is important to know how the coefficients behave under some linear transformations in the original galaxy image, 
such as rotations, dilation or contractions or an increasing of the fiux. due to the linearity of the CHEF functions, 
calculation the translation of these transformations into CHEF domain is straightforward. 

• Rotation. Given f^°*{r,0) :— f{r,6 + p) with p e [— 7r,7r] any angle, it can be proved that 

rrot r ci^^l-P 

J nm ~ J nm e 

That is, a rotation of p radians in real space is equivalent to a rotation of mp radians in CHEF space. 

• Dilation or contraction. Given a £ M and defining f^^^rjO) :— f{ar,9), it yields 

j-dil J- 

Jnm;L ~ Jnm;aL- 

What indicates that a dilation or contraction by a factor a in real space is equivalent to a dilation or contraction 
by the same factor in the parameter scale L in CHEF space. 

• Flux increasing. Given fc G M and /™'^(r', 6) := k ■ f{r, 9), it is obtained 

fine _ I, £ 
J nm ~ Jnm- 

This means that an increasing in the flux by a factor k is equivalent to an increasing in the CHEF coefficients 
by the same factor. 
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